from util import *
from glob import glob
/usr/local/lib/python3.8/dist-packages/geopandas/_compat.py:124: UserWarning: The Shapely GEOS version (3.11.2-CAPI-1.17.2) is incompatible with the GEOS version PyGEOS was compiled with (3.9.1-CAPI-1.14.2). Conversions between both will be slow. warnings.warn(
df = load_AOIs()
df
| Taranaki | AOI | SSP 4.5 (p50) | SSP 4.5 (p83) | SSP 8.5 (p50) | SSP 8.5 (p83) | Rate SSP 4.5 (p50) | Rate SSP 4.5 (p83) | Rate SSP 8.5 (p50) | Rate SSP 8.5 (p83) | match | match_score | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 7 | NORTH | TongaporutuRiver | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | TongapurutuRiverCliffs | 93.750000 |
| 11 | SOUTH | HangatahuaRiver_South | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | HangatahuRiver_South | 97.435897 |
| 21 | SOUTH | Rahotu | 0.58 | 0.78 | 0.84 | 1.10 | 0.0058 | 0.0078 | 0.0084 | 0.0110 | Rahotu | 100.000000 |
| 20 | SOUTH | Pihama | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | Pihama | 100.000000 |
| 19 | SOUTH | OpunakeBeach | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | OpunakeBeachCliffs | 100.000000 |
| 18 | SOUTH | OhaweBeach | 0.57 | 0.78 | 0.83 | 1.10 | 0.0057 | 0.0078 | 0.0083 | 0.0110 | OhaweBeach | 100.000000 |
| 17 | SOUTH | Oeo | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | Oeo | 100.000000 |
| 16 | SOUTH | Manutahi | 0.57 | 0.78 | 0.83 | 1.10 | 0.0057 | 0.0078 | 0.0083 | 0.0110 | Manutahi | 100.000000 |
| 15 | SOUTH | ManaBay | 0.57 | 0.78 | 0.83 | 1.10 | 0.0057 | 0.0078 | 0.0083 | 0.0110 | ManaBayCliffs | 100.000000 |
| 14 | SOUTH | KaupokonuiBeach | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | KaupokonuiBeach | 100.000000 |
| 13 | SOUTH | Kakaramea | 0.57 | 0.78 | 0.83 | 1.10 | 0.0057 | 0.0078 | 0.0083 | 0.0110 | Kakaramea | 100.000000 |
| 12 | SOUTH | Hawera_WaihiBeach | 0.57 | 0.78 | 0.83 | 1.10 | 0.0057 | 0.0078 | 0.0083 | 0.0110 | Hawera_WaihiBeach | 100.000000 |
| 0 | NORTH | Mohakatino_River | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | MohakatinoRiver | 100.000000 |
| 10 | SOUTH | CapeEgmont | 0.58 | 0.78 | 0.84 | 1.10 | 0.0058 | 0.0078 | 0.0084 | 0.0110 | CapeEgmont | 100.000000 |
| 9 | NORTH | Waitara | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | Waitara | 100.000000 |
| 8 | NORTH | UrenuiRiver_North | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | UrenuiRiverNorthCliffs | 100.000000 |
| 6 | NORTH | PariokariwaPoint | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | PariokariwaPointCliffs | 100.000000 |
| 5 | NORTH | Onaero | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | OnaeroCliff | 100.000000 |
| 4 | NORTH | Oakura_South | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | OakuraSouth | 100.000000 |
| 3 | NORTH | Oakura | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | Oakura | 100.000000 |
| 2 | NORTH | NewPlymouth_Airport | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | NewPlymouthAirport | 100.000000 |
| 1 | NORTH | NewPlymouth | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | NewPlymouthCliffs | 100.000000 |
| 22 | SOUTH | WainuiBeach | 0.57 | 0.78 | 0.83 | 1.09 | 0.0057 | 0.0078 | 0.0083 | 0.0109 | WainuiBeach | 100.000000 |
| 23 | SOUTH | WaipipiBeach | 0.57 | 0.78 | 0.83 | 1.09 | 0.0057 | 0.0078 | 0.0083 | 0.0109 | WaipipiBeachCliffs | 100.000000 |
site = df.match.sample(1).iloc[0]
site
'Oeo'
gdf = gpd.read_file(f"Shapefiles/{site}_intersects.shp")
gdf.crs = 2193
gdf = enrich_df(gdf)
transect_metadata = get_transect_metadata(f"Shapefiles/{site}_TransectLines.shp")
if site == "ManaBayCliffs":
print("Flipping")
for k, v in transect_metadata.items():
transect_metadata[k]["Azimuth"] = v["Azimuth"] + 180
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance return lib.distance(a, b, **kwargs)
linear_models = fit(gdf, transect_metadata)
linear_models.loc[linear_models.slope > 0, "slope"] = pd.NA
rolled_slopes = linear_models.groupby("group").slope.rolling(10, min_periods=1).mean().dropna().reset_index(level=0)
linear_models.slope = rolled_slopes.slope
linear_models.dropna(inplace=True)
linear_models
| TransectID | slope | intercept | group | |
|---|---|---|---|---|
| 0 | 2 | -0.025895 | 11.839523 | 0 |
| 1 | 3 | -0.026103 | 11.394666 | 0 |
| 2 | 4 | -0.032111 | 12.067457 | 0 |
| 3 | 5 | -0.037306 | 12.952147 | 0 |
| 4 | 6 | -0.041094 | 14.552892 | 0 |
| ... | ... | ... | ... | ... |
| 665 | 806 | -0.088458 | 17.144344 | 50 |
| 666 | 807 | -0.098385 | 17.888823 | 50 |
| 667 | 808 | -0.100127 | 16.604210 | 50 |
| 668 | 809 | -0.105771 | 5.736368 | 50 |
| 669 | 810 | -0.105771 | 4.466733 | 50 |
657 rows × 4 columns
results = predict(gdf, linear_models, transect_metadata)
results
| TransectID | BaselineID | group | Year | ocean_point | linear_model_point | sqrt_model_point | BH_model_point | Sunamura_model_point | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 2.0 | 1 | 0.0 | 2100 | POINT (1687908.5767938623 5620307.34597674) | POINT (1688034.2711515864 5620793.369995707) | POINT (1688035.5534006841 5620798.3280850835) | POINT (1688043.1672309272 5620827.768582187) | POINT (1688037.0176938204 5620803.9900865285) |
| 1 | 3.0 | 1 | 0.0 | 2100 | POINT (1687915.1713428975 5620306.442715614) | POINT (1688044.0861092217 5620791.621524045) | POINT (1688045.4118067573 5620796.610869867) | POINT (1688053.2101313008 5620825.9603533875) | POINT (1688046.9030184844 5620802.223138728) |
| 2 | 4.0 | 1 | 0.0 | 2100 | POINT (1687923.1177022296 5620304.70048816) | POINT (1688053.6788356411 5620788.917562359) | POINT (1688055.3321562614 5620795.049295843) | POINT (1688062.928654537 5620823.222719841) | POINT (1688056.5343737882 5620799.508006012) |
| 3 | 5.0 | 1 | 0.0 | 2100 | POINT (1687929.0214491673 5620302.536182147) | POINT (1688063.1581929883 5620785.8900468) | POINT (1688065.1311122181 5620792.999345171) | POINT (1688072.6592310544 5620820.126478494) | POINT (1688066.0911000874 5620796.458604978) |
| 4 | 6.0 | 1 | 0.0 | 2100 | POINT (1687943.6634779107 5620297.272404213) | POINT (1688072.646862405 5620782.368306764) | POINT (1688074.7352682156 5620790.222629285) | POINT (1688081.7768776012 5620816.705543141) | POINT (1688075.4651071057 5620792.967493336) |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 652 | 806.0 | 1 | 50.0 | 2100 | POINT (1687794.0844616203 5620310.788355438) | POINT (1687889.1002257578 5620803.734863906) | POINT (1687892.4113368606 5620820.913072934) | POINT (1687895.824934831 5620838.622988925) | POINT (1687891.1748046514 5620814.497882353) |
| 653 | 807.0 | 1 | 50.0 | 2100 | POINT (1687821.543196849 5620306.082154956) | POINT (1687899.1019846308 5620802.495838878) | POINT (1687902.105572404 5620821.720250634) | POINT (1687904.5866393007 5620837.600276526) | POINT (1687900.7938008038 5620813.324279123) |
| 654 | 808.0 | 1 | 50.0 | 2100 | POINT (1687873.6267297296 5620298.320708633) | POINT (1687909.2520916292 5620802.726387996) | POINT (1687910.64721948 5620822.479462039) | POINT (1687911.7553046227 5620838.168409724) | POINT (1687910.0242252424 5620813.6587283) |
| 655 | 809.0 | 1 | 50.0 | 2100 | POINT (1687893.8018985782 5620296.623321697) | POINT (1687919.8082378658 5620813.567189055) | POINT (1687920.859272824 5620834.459251945) | POINT (1687921.5934341024 5620849.052623445) | POINT (1687920.3588562713 5620824.512166921) |
| 656 | 810.0 | 1 | 50.0 | 2100 | POINT (1687939.3591070531 5620298.028960121) | POINT (1687928.690593155 5620815.9733632095) | POINT (1687928.2598098365 5620836.8874109285) | POINT (1687927.9589022126 5620851.496139112) | POINT (1687928.4649134837 5620826.929858534) |
657 rows × 9 columns
gpd.GeoSeries(results.linear_model_point).distance(gpd.GeoSeries(results.Sunamura_model_point)).describe()
count 657.000000 mean 10.961338 std 0.009703 min 10.924303 25% 10.957830 50% 10.964826 75% 10.967859 max 10.972655 dtype: float64
results.set_geometry(f"linear_model_point", inplace=True, crs=2193)
poly = prediction_results_to_polygon(results)
poly.explore(tiles="Esri.WorldImagery")
site = df.match.sample(1).iloc[0]
site = "CapeEgmont"
print(site)
gdf = gpd.read_file(f"Shapefiles/{site}_intersects.shp")
gdf.crs = 2193
gdf = enrich_df(gdf)
transect_metadata = get_transect_metadata(f"Shapefiles/{site}_TransectLines.shp")
m = gpd.read_file(f"Shapefiles/{site}_TransectLines.shp").explore("TransectID", tiles="Esri.WorldImagery", max_zoom=22)
linear_models = fit(gdf, transect_metadata)
linear_models.loc[linear_models.slope > 0, "slope"] = pd.NA
rolled_slopes = linear_models.groupby("group").slope.rolling(10, min_periods=1).mean().dropna().reset_index(level=0)
linear_models.slope = rolled_slopes.slope
linear_models.dropna(inplace=True)
results = predict(gdf, linear_models, transect_metadata)
results["geometry"] = results["linear_model_point"]
pd.concat([gdf, results])[["Year", "geometry"]].explore("Year", m=m)
gpd.GeoSeries(results.ocean_point, crs=2193).explore(m=m, color="blue")
#gpd.GeoSeries(poly, crs=2193).explore(color="pink", m=m)
results.set_geometry(f"linear_model_point", inplace=True, crs=2193)
polygon = prediction_results_to_polygon(results)
polygon.explore(m=m)
m
CapeEgmont
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance
return lib.distance(a, b, **kwargs)
<ipython-input-10-8352519b00ac>:15: FutureWarning: You are adding a column named 'geometry' to a GeoDataFrame constructed without an active geometry column. Currently, this automatically sets the active geometry column to 'geometry' but in the future that will no longer happen. Instead, either provide geometry to the GeoDataFrame constructor (GeoDataFrame(... geometry=GeoSeries()) or use `set_geometry('geometry')` to explicitly set the active geometry column.
results["geometry"] = results["linear_model_point"]
run_all_parallel()
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance return lib.distance(a, b, **kwargs) /usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance return lib.distance(a, b, **kwargs) /usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance return lib.distance(a, b, **kwargs) /usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance return lib.distance(a, b, **kwargs) /usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance return lib.distance(a, b, **kwargs) /usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance return lib.distance(a, b, **kwargs) /usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance return lib.distance(a, b, **kwargs) /usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance return lib.distance(a, b, **kwargs) /usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance return lib.distance(a, b, **kwargs) /usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance return lib.distance(a, b, **kwargs) /usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance return lib.distance(a, b, **kwargs) /usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance return lib.distance(a, b, **kwargs) /usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance return lib.distance(a, b, **kwargs) /usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance return lib.distance(a, b, **kwargs) /usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance return lib.distance(a, b, **kwargs) /usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance return lib.distance(a, b, **kwargs) /usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance return lib.distance(a, b, **kwargs) /usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance return lib.distance(a, b, **kwargs) /usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance return lib.distance(a, b, **kwargs) /usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance return lib.distance(a, b, **kwargs)
Flipping
/usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance return lib.distance(a, b, **kwargs) /usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance return lib.distance(a, b, **kwargs) /usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance return lib.distance(a, b, **kwargs) /usr/local/lib/python3.8/dist-packages/pygeos/measurement.py:70: RuntimeWarning: overflow encountered in distance return lib.distance(a, b, **kwargs)
def read_file(f):
df = gpd.read_file(f)
df["filename"] = f
return df
samples = pd.concat(read_file(f) for f in glob("Projected_Shoreline_Polygons/*.shp"))
samples.explore("filename", tiles="Esri.WorldImagery", style_kwds={"fill": False})